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ABSTRACT 


The power spectrum for a stationary random process can be 
defined with the Wiener-Khintchine Theorem, which says that the power 
spectrum and the autocorrelation function are a Fourier transform 
pair. To implement this theorem for signals that are discrete and of 
finite length we can use the Blackman-Tuckey method. Blackman and 
Tukey (1958) show that a function w(r), called a lag window, can be 
applied to the autocorrelation estimates to obtain power spectrum 
estimates that are statistically stable. The Fourier transform of w(r) 
is called a spectral window. 

Typical choices for spectral windows show a distinct trade-off 
between the mainlobe width and sidelobe strength. A new idea for 
designing windows by taking linear combinations of the standard 
windows to produce hybrid windows was introduced by Smith (1985). We 
implement Smith's idea to obtain spectral windows with narrow 
mainlobes and smaller (compared with typical windows) near sidelobes. 

One of the main contributions of this thesis is that we show 
that Smith's problem is equivalent to a Quadratic Programming (QP) 
problem with linear equality and inequality constraints. A computer 
program was written to produce hybrid windows by setting up and 
solving the QP problem. We also developed and solved two variations of 
the original problem. The two variations involved changing the 
inequality constraints in both cases from nonnegativity on the 
combination coefficients to nonnegativity on the hybrid lag window 
itself. For the second variation, the window functions used 
to construct the hybrid window were changed to a frequency-variable 


ix 


set of truncated cosinusoids. 

A series of tests was run with the three computer programs to 
investigate the behavior of the hybrid spectral and lag windows. 
Emphasis was put on obtaining spectral windows with both relatively 
narrow mainlobes and the lowest possible (for these algorithms) near 
sidelobes. Some success was achieved for this goal. A 10 dB peak 
sidelobe reduction over the rectangular spectral window without 
significant mainlobe broadening was achieved. Also, average sidelobe 
levels of -117 dB were reached at a cost of doubling the mainlobe 
width (at the -3 dB point). 


X 


1 


tri 


CHAPTER 1 
INTRODUCTION 

The definition of the power spectrum P(f) due to Wiener 
(Robinson, 1980) says that P(f) forms a Fourier transform pair with a 
quantity called the autocorrelation function. The autocorrelation of a 
stationary, ergodic realization of a random process (Blackman and 
Tukey, 1958) is 


C(r ) - Lim 1/T 
T -► 00 


f T/2 

x(t)x(t+r) dt. 

J -t/2 


T is the length of a segment of x(t) and r is a variable called the 
lag time. It corresponds to the amount of time shift between x(t) and 
x(t+r) when computing the auto (self) correlation. The power spectrum 
P(f) is the Fourier transform of C(r): 


P(f) - 


-i2*fr 

C(r)e dr 


where (f) is the frequency. To implement this definition of P(f) for a 
finite interval of a sampled time series, we need discrete estimators 
of C(r ) and the Fourier transform (Oppenheim and Schafer, 1975). For 
a real (N) point sequence x[n] , we can write the autocorrelation 


N- 1 m | 

c[m] - 1/(N- |m| ) 


-1 

x[n]x[n+m] . 


n— 0 


function as 



2 


and 


( 1 . 1 ) 



■ i2wkn/N 


m—0 


for the discrete Fourier transform (DFT) of the autocorrelation (power 
spectrum). Blackman and Tukey (1958) realized that the application of 
the above two formulae leads to statistically unstable power spectrum 
estimates. The variance of the estimate can be quite large and does 
not decrease with increasing (N) . Blackman and Tukey solved this 
problem by noticing that the later lags (m) of c[m] have a smaller 
number of products to average over, and hence, are less reliable 
estimates. They showed that spectrum estimates made by ignoring later 
lags became more stable. For Gaussian processes, they recommend 
keeping only about 10% of the autocorrelation lags. The theoretical 
properties of lag windows are more conveniently discussed for the 
continuous case so we will stay with the variable t for the rest of 
this chapter. 

Deleting later lags amounts to multiplying the autocorrelation 
estimate with a rectangle function Il(r/L) of the appropriate length 
(L) , where 


n(r/L) 


f 1, M < L/2 


l 0, | r | > L/2 . 

Note that (Bracewell, 1978) II(r/L) D L sinc(Lf) - sin(irLf )/(irf ) (where 
D denotes the Fourier Transform) . According to the Convolution Theorem 
(Bracewell, 1978), this function, sin(irLf)/(irf) , will then be 


convolved with our spectral estimate. Figure 1 shows this function for 


L - 1 second. Note that the peak sidelobe is at about -13.5 dB. Our dB 


3 


convention uses the following definition 

W(f) (dB) - 20 log 1Q | w ( f )/ w max l • 

The convolution of the rectangle transform with our spectrum will 
reduce the variance of our estimate but will introduce other problems 
at the same time. First, note the main lobe width. This will 
fundamentally limit our resolution in the spectrum. We are especially 
concerned about this if we are trying to resolve closely spaced peaks . 
Also, the sidelobes themselves can mask weak signals in the presence 
of strong resonances, and generally distort the spectrum (Marple, 

1987). Side lobe distortion is sometimes referred to as leakage. 

Hence, we need to pay close attention to these effects. 

First, we make the observation that functions other than II(r/L) 
should also serve to reduce variance (Blackman and Tukey, 1958). 
Consider a function w(r) , called a lag window, which has the following 
properties . 

1) w(r ) > 0, for all r 

2) w(r ) - 0 , | r | > L/2 

3) w(r ) piecewise continuous 

4) w(-r) — w(r) (symmetric about the origin) 

5) w(0) - 1 
Note that 

, T/2 

C(0) - Lim 1/T [x(t)] dt. 

T -* ® I 

J -T/2 

C(0) is called the power of the signal. Property 5) simply 
insures that we preserve the power. Note also that the inverse 


4 


mu 


. fca 


transform is 


i2wfr „ 
C(r ) - P(f)e df, 


so 


C(0) 


00 

P(f) df. 

-00 


Thus, property 5) preserves the area under the power spectrum P(f) . 
Property 4) guarantees that P' (f) , the Fourier transform of the 
windowed autocorrelation, is real, where 


P'(f) - 


C(r )w(r)e 


■ i2*rfr 


dr , 


Let w(r ) 3 W(f) . 

We call W(f) the spectral window. If w(r) is real and even, W(f) is 
real and even, since 


W(f) 


r 00 

w(r )cos(2*fr ) 
-00 


dr + i 




<*> 

w(r)sin(2*rr ) 
00 


dr. 


An even function w(r) multiplied by the odd function sin(2?rfr) is odd, 
and the integral over r vanishes. Now, P'(f) - P(f)*W(f) . The 
convolution of two real and even functions is real and even (P(f) is 
real and even for real x(t)>. 

Eleven typical spectral windows, common (except for window 4) in 
the literature, are shown in Figures 1-22 (Smith, 1985; Kreamer, 1988; 
Marple , 1987). Refer to Table 1 for the formulae. Each window has the 
following properties which are useful to summarize overall behavior 



(Marple, 1987); 

1) ma inlobe width (typically referenced to the -3 dB point), 

2) sidelobe decay rate (dB/octave) , 


5 


and 

3) the peak sidelobe (amplitude of the largest sidelobe) . 

When observing spectral window behavior in Figures 1-22, it is 
immediately obvious that there is a trade-off between mainlobe width 
and the peak sidelobe. Typically, having small (near) sidelobes and 
narrow mainlobes is a conflicting goal (Marple, 1987). One problem in 
spectral analysis where this trade-off is particularly bothersome is 
dicussed next. 

Consider a signal x(t) that contains a large amplitude 
sinusoidal component of frequency f ' . This "resonance" would ideally 
show up in the power spectrum as a peak at frequency f ' . However, our 
spectral window W(f) will be convolved with the peak, and the 
sidelobes may mask weak signals at frequencies near f ' . Our spectral 
window must have a mainlobe narrow enough to resolve the two signals, 
and needs small near sidelobes if we are to detect the weak spectral 
response. This creates a challenging window design problem in trying 
to overcome these apparently conflicting goals. The problem of 
designing such a spectral window will constitute the main emphasis of 
the rest of this work. 

Smith (1985) posed a mathematical statement of a possible scheme 
to design the desired window. He suggested that we might be able to 
generate a better window than the standard literature windows if we 
take a linear combination of them. Let's define the windows in Table 1 


6 


as w^r), where i-l,M (1 to 11 here). The linear combination will 

produce a "hybrid" window w (r). Then, 

H 


M 


V° ■ / Vi (r) - 

i-1 


To insure that w (r) has property 1), we can impose a nonnegativity 
H 


condition on the a^, giving 


a i 0 . 
i 


To insure property 5) we need 


M 


•i- 1 - 


Note that if 
then, since 
we have 
where 


i-1 


w(r) D V f) ' 

H n 

w (r) - 2 a i w i (r), 
W R (f) - 2 a^Cf), 
w (r) 3 W i (f). 


For a criterion to determine an appropriate set of a^, we can use the 
least squares measure to insure that the spectral window is small over 
a desired frequency band f^ to f^. That is, choose a^ such that 


.f 


I W u (f) | df 


H 


is a minimum. This problem, and two variations to be discussed later, 
form the basic subject of study in this work. A brief survey of the 
contents of the rest of this thesis follows next. 
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In Chapter 2 we proceed to solve the problem of how to determine 

the a so that w is a valid lag window while minimizing the above 
i H 

integral. We will also discuss two variations of the original problem. 
In Chapter 3 we discuss the algorithmic details of implementing (as 
computer codes) the solutions presented in Chapter 2, and discuss a 
series of tests that were run with the programs. The figures (1 
through 80) that follow after Chapter 4 are also discussed. In Chapter 
4 some general conclusions are drawn about algorithmic behavior from 
the tests that we discussed in Chapter 3. 


8 


CHAPTER 2 
THEORY 

In Chapter 1 we posed an optimization problem with the following 


form 




min F(a) - 


2 2 
I W (f) | df, 

ri 


where 


subject to a^ > 0 
M 


and 


‘i- 1 - 


i-i 


M 


W (f) - 
H 


a W.(f) 
i l 


i-1 


and 


Since 


* - Ia i V- 

w D W t (f) 


and the w are even, we have that the W^f) are real. Thus, 


F(a) 


- C 

-f 


(f) 1 df 


2 2 ^ 
[| a i W t (f)] df 

1 


- ? ? h a j 


H (f) W (f) df 

£ 1 J 
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Now we know the W^(f) so, in principle, we could carry out the 

indicated integration. In practice, this would require numerical 

quadrature techniques. Another approach would be to evaluate W (f) 

i 

numerically with the discrete Fourier transform (DFT) , and change the 
integral to a discrete sum over W^(k)W (k). Recall that the DFT was 
defined in Chapter 1 in Eq. No. (1.1). We will proceed in the following 
manner. Let 




W.(f) W^(f) df 


W^k) W^k) 


k-k 


or 


W W 
ki kj 


k-k 


where k^ and k^ are the indices that correspond to the frequency 


points f^ and f . That is, 


f - (k - 1) Af, 
n n 

where Af - l/(NAr), and A r is the sample rate for C(r) (to give c[m]). 
N is the total number of sample points. 

Consider a matrix W with elements W , where (i) is the window 

ij J 

number and (i) is the frequency index, tf is (k^ - k + 1) X M in size. 
Then , 


2 


w w 

ki kj 





t l th 

where W,, is an element of W and S W W is the ii element of 
ik k ki kj J 

T T T 

W W. W W is square, M X M, and symmetric. Let H - ¥ V for convenience. 


' ' : rni nr I % 

i Nlliii J 
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El 


U 


i t3 


Then 

F(.) - | J V . | « ki » kj - l J Vl/j 

-f Vf 

which can be written as 

F(a) - a T Ha - a T (V T V)a - (Wa) T (Wa) 

In general, a quadratic form (Scales, 1985) is a function 
of x that can be written as 

T T 

f(x) - 1/2 x Ax + b x + c. 

A is called the Hessian. Obviously, our objective function F(a) is a 
quadratic form, where we identify 

1/2 A - H 
b - 0 
c - 0. 


We have immediately that H is symmetric, since 


or 


T 

H - W V 

T T T T 

H - (W V) - V W - H. 

T 

H - H . 


Furthermore, F(a) 


2 2 
| W (f) | df i 0, 

c H 


so F(a) > 0 for any choice of a. By definition, a matrix A is 

T 

nonnegative definite if and only if, for any x, x Ax i 0. Hence, since 

T 

F(a) - a Ha > 0 

for any a, we have that H is symmetric and nonnegative semidef inite . 

It is instructive to note that the vector a that minimizes F(a) is 
(0,...,0), or the trivial solution. This is shown in Appendix A. 


,in li b 
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In summary, then, we wish to find an a that miminizes 

T 

a Ha 

subject to ^ a^ - 1 

and a, >0, 

l 

given that H is symmetric and nonnegative semidef inite . This 

is a special case of a more general mathematical problem in 

optimization theory call the Quadratic Programming problem (Bronson, 

1982) . The full formulation allows more general linear constraints and 

a nonzero b and c (although c does not affect the solution) . 

Quadratic programming (QP) is an effective way In general to 

deal with least squares problems if inequality constraints are needed. 

Various approaches for the solution of this problem are discussed by 

Bronson (1982), Lawson and Hanson (1974), Scales (1985), and Hillier 

and Lieberman (1967). The solution to the QP problem investigated In 

T M 

this work is a code, called DQPROG, from the IMSL MATH /Library. The 
details of QP are beyond the scope of the work, but a brief descrip- 
tion is given in Appendix B. 

In addition to the QP problem we have just outlined, two 
variations were also studied in this work. We will refer to 
the above QP problem henceforth as Program- 1. The first variation we 
will consider Is a relaxation of the nonnegativity constraints on a. 
Instead, we will impose this condition on the hybrid lag window 
itself. The inequality constraints then have the following form: 

Ba > 0, 

th th 

where B Is the lag window value for the i time sample and the j 

window. B is (N X M) , where N is the number of lag window samples used 
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for the DFT. This new problem will be referred to as Program- 2. 

The second variation on Program- 1 we wish to consider is that of 

retaining the constraints of Program- 2, but changing the functions 

w (r) . The new functions are (Nuttall, 1981) 
i 

(2.1) w - cos(2?rLr/T) , L - 0 , 1 , . . . 

L+l 

where T is the window length and L is an integer . They offer the 
advantage of being more general than those of Program-1, and are 
reminiscent of a Fourier construction. This variation will be called 
Program- 3. Examples of the implementation of Programs-1, 2 and 3, with 
a discussion of results, are presented in Chapter 3. 
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CHAPTER 3 
RESULTS 

Introduction 

A general outline of the algorithm for Program- 1 is given first, 
followed by an introduction to the tests that were run. Program- 2 and 
Program- 3 follow this outline as well; only the constraints and/or the 
window functions differ. 

1) The Discrete Fourier Transform (DFT) was discussed in Chapter One. 
An algorithm for efficient evaluation of the DFT is called the Fast 
Fourier Transform (FFT) . An FFT (Claerbout, 1976) of the radix 2 type 
was used to evaluate W. The w^(r) (see Table 1 and Figures 1-22) were 
sampled at 256 points, and zero padded to 1024 points before computing 
the FFT. The zero padding was used to obtain a more densely sampled 
transform. 

2) W was constructed by using the window FFT's over the chosen 

frequency band. The FFT over this band for window w^(r) becomes the 

th T 

1 column of W. Then H was constructed by H - O. 

H and the constraints were given to DQPROG (an IMSL quadratic 
programming subroutine), which returned a solution vector a. 

3) The hybrid spectral window was constructed in the function domain 
using ^ a^w^(r) , and then sampled and the FFT computed in the same 
manner as in step 1 (except that the zero padding was to 2048 points) . 

Now we will consider the testing that was done using Programs-1, 

2 and 3. These three programs were tested for performance by varying, 
in each case, the frequency band we wish to attenuate, computing two 
quantitative measures that relate to performance (discussed below) , 
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and displaying the resulting hybrid spectral windows along with a 
graph of the rectangular spectral window for comparison. The two 
quantitative measures of performance are the mean dB level (MDB) over 
the desired attenuation band and the half power width (HPW) . The half 
power width is the frequency (Hz) that corresponds to the -3 dB point 
on the spectral window of interest. The mean dB levels and half power 
widths are listed for each test in Table 2 through Table 5. 

The testing strategy was as follows. For Program- 1 and 
Program- 2, the eleven lag windows listed in Table 1 and Figures 1 
through 22 were used to compute the hybrid spectral windows. For 
Program- 3, it was of interest to see how the results would vary with 
different numbers of the cosine windows, since the window number 
scheme (window number 1,2, etc.) corresponds to the frequency of the 
lag window. That is, increasing window numbers mean increasing 
frequencies (see Figures 23 through 40). The maximum number of cosine 
windows tested was 9 (corresponding to L - 8 in equation number 2.1). 

For all three programs, the desired frequency band for 
attenuation was varied. Figures 41 through 50 show results for the 
frequency band 1 through 2 Hz. Figures 51 through 60 show results for 
the frequency band 2 through 4 Hz. Figures 61 through 70 show results 
for the frequency band 2 through 8 Hz. Figures 71 through 80 show 
results for the frequency band 4 through 8 Hz. The rationale for this 
testing scheme was to observe window behavior for the cases in which 
the frequency band was close to and farther away from the origin, and 
in which the frequency band width itself varied. 

In the case of Program- 3, for each frequency band tested, three 
different combinations of windows were used. The combinations were. 
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window numbers 1 through 3 (L — 0 through L ™ 2) , window numbers 1 
through 6 (L - 0 through L - 5), and window numbers 1 through 9 (L - 0 
through L - 8) . 

Discussion of Tests 
Test for 1 through Z Hz 

For the 1 - 2 Hz test (Figures 41 through 50), the MDB's range 
from -27 dB to -29 dB, and the HPW's range from 0.46 to 0.48 Hz (refer 
to Table 2). The average MDB is about -28 Hz and the average HPW is 
about 0.47 Hz. These values are fairly clustered, so to further judge 
performance, we can visually inspect the spectral and lag windows in 
Figures 41 through 50. For comparison purposes, the rectangular spec- 
tral window (also shown in Figure 1) is overplotted with dashed lines. 
Note that the HPW of the rectangular spectral window (RSW) is 0.44 Hz, 
and the peak sidelobe is at -13.5 dB. 

For each case, we note that the HPW's of the hybrid spectral 
windows are relatively close to the RSW width, while the peak 
sidelobes of the hybrid windows inside the attenuation band are at 
about -23 dB. This is a 10 dB improvement over the RSW without much 
gain in window width (increases by a factor of 1.07), which is an 
encouraging result. 

The relative performance of Programs 1, 2, and 3 (3 windows, 6 
windows, and 9 windows) is discussed next. First, let's adopt a short- 
hand notation. Let Program- 1 be represented by PI, Program- 2 by P2 , 
and Program- 3, 3 windows, 6 windows, and 9 windows be respectively 
represented by P3 ; 3 , P3 ; 6 , and P3 ; 9 . We can summarize the relative 
window behavior, which we observe visually, by saying that the side- 
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lobe behavior, outside of the attenuation band, of P2 (Figure 43), 

P3 ; 6 (Figure 47), and P3;9 (Figure 49) is relatively poor. On the 
other hand, PI (Figure 41) and P3;3 (Figure 45) have sidelobe behavior 
outside of the attenuation band that is quite near that of the RSW. As 
we would expect, the lag windows for PI and P3;3 (Figures 42 and 46) 
look quite similar, while the lag windows for P2, P3;6, and P3;9 
(Figures 44, 48, and 50) are quite variable with respect to each 
other. It is interesting to compare the lag windows for PI and P3;3 
with the rectangular lag window (Figure 2) . The hybrid lag windows 
start at 1.0, smoothly taper off to about 0.6, and then stay fairly 
level. This behavior is somewhat different from the standard lag 
windows (included in Figures 2 through 22). 

Test for 2 through 4 Hz 

For the 2 - 4 Hz test (Figures 51 through 60), the MDB's range 

from -54 dB to -64 dB, and the HPW's range from 0.62 through 0.66 Hz 

(refer to Table 3). The average MDB is about -58 dB and the average 
HPW is about 0.64 Hz. Note that there is more variability in the MDB's 
than was the case for the 1 - 2 Hz test. However, the HPW's do not 
significantly vary. 

In moving the frequency band away from the origin (0 Hz) , we 
have allowed the average HPW to increase over the RSW width (0.44 Hz). 
The Increase Is by a factor of 1.45. Also, by moving the attenuation 
band, we have improved our attenuation to an average dB level of -58 
dB. 

In comparing PI (Figure 51) to P2 (Figure 53) , we note that P2 

has a slightly better attenuation level in the 2 - 4 Hz band, but 

larger side lobes outside of the band. However, P2's behavior outside 
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of the attenuation band is better than in was in the 1 - 2 Hz case. PI 
and P3;3 (Figure 54) have very similar behavior. P2 has slightly 
better behavior inside the attenuation band than P3;3, and slightly 
worse behavior outside of the band. P3;6 (Figure 57) and P3;9 (Figure 
59) exhibit similar behavior in the 2 - 4 Hz band, and have comparable 
mainlobe widths, but P3;6 has much better sidelobe behavior outside of 
the attenuation band than does P3;9. Compared to the other windows 
(except P3;9) however, P3;6 has poor sidelobe behavior outside of the 
2 - 4 Hz band. Note that the lag windows for PI (Figure 52), P2 
(Figure 54), and P3 ; 3 (Figure 56) are quite similar. The lag windows 
for P3;6 (Figure 58) and P3;9 (Figure 60) are considerably different 
from each other and the other lag windows. The P3;9 lag window has an 
oscillatory behavior uncharacteristic of typical lag windows. However, 
P3;9 has the best attenuation characteristics for the 2 - 4 Hz case 
in the attenuation band. 

Test for 2 through 8 H£ 

For the 2 - 8 Hz test (Figures 61 through 70), the MDB's range 
from -55 dB to -57 dB, and the HPW's range from 0.62 Hz to 0.66 Hz 
(refer to Table 4). The average MDB is about -56 dB and the average 
HPW is about 0.67 Hz. It is evident from Figures 62 through 69 that 
the results for the spectral and lag windows are fairly uniform as 
compared with the last two tests. Going to a wider bandwith seems to 
have a stabilizing effect on the programs' behavior. It would be 
difficult to choose one window over the other from this test, except 
for the P3 ; 9 (Figure 69) window, which has the poorest behavior and 
would be excluded. Towards the end of the frequency band, the 
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sidelobes of P3;9 come back to a higher level than for the other 
spectral windows. 

In comparison to the 2 - 4 Hz case, the Pi's (Figures 51 and 
61) are about the same. For the P2's (Figures 53 and 63), a 
significant improvement in the behavior of the 2 - 8 Hz case can be 
noted. The sidelobes from 4 - 8 Hz are much lower than for the 2-4 
Hz case. This makes a lot of sense, and shows that the algorithms are, 
fortunately, behaving in an intuitive and reasonable way. The results 
are similar but even more dramatic when comparing the P3;6 and P3;9 
cases for the 2 - 8 Hz test (Figures 67 and 69) with the 2 - 4 Hz test 
(Figures 57 and 59) . It seems safe to conclude that we gain more than 
we lose by opening up the bandwidth when possible. 

Test for 4 through 8 Hi 

For the 4 - 8 Hz test (Figures 71 through 80), the MDB's range 
from -83 dB to -117 dB, and the HPW's range from 0.81 Hz to 0.93 Hz 
(refer to Table 5). The average MDB is about -101 dB, and the mean HPW 
is about 0.88 Hz. The average MDB has almost doubled over the 2 - 8 Hz 
case, and the average HPW has increased by a factor of about 1.3. We 
note for this test that we have a broader range of results in terms of 
attenuation levels and mainlobe widths across the various programs. 

Specifically, P3;6 (Figure 77) and P3;9 seem to be in a class by 
themselves with much higher attenuation levels (better than -100 dB) , 
but broader main lobes. However, P3;6 and P3;9 are fairly comparable 
to each other. It is also interesting to note that sidelobe behavior 
outside the attenuation band is fairly good for all cases. The 4-8 
Hz band, being farther from the origin, seems to allow more stable 
results. This is also reflected in the lag window behavior. The lag 
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windows for the 4 - 8 Hz case (Figures 72, 74, 76, 78, and 80) are all 
fairly comparable and well behaved (smoothly tapering to a low value) . 

The spectral windows for PI (Figure 71) , P2 (Figure 73) , and 
P3;3 (Figure 75) are similar. Their behavior can be characterized in 
the following way. As we progress through the above sequence, the 
overall mean attenuation improves, but the attenuation nearest the 
mainlobe deteriorates. We can summarize the extremes of this test by 
saying that PI has the narrowest mainlobe behavior, and P3;9 has the 
best sidelobe behavior. One would have to choose from them according 
to need. 

Conclusions 

We will summarize here some of the conclusions that are 
suggested by the above tests. The behavior differences between 
Programs- 1, 2, and 3 can be significant. This is illustrated in the 1 
- 2 Hz test by the instabilities (large sidelobes) exhibited for 
Program- 3, and the differences between Program- 1 and Program- 3 results 
in the 4 - 8 Hz test. Program-1 and Program-2 show their greatest 
differences for cases where the attenuation band is close to the 
origin. 

In general, Program- 2 seems to have less stable sidelobe 
behavior outside the attenuation band than does Program-1. In most of 
the cases that we looked at, Program- 1 and Program- 3 (3 windows) 
seemed to be the most similar, but for the case away from the origin 
(4-8 Hz), Program-3 (6 windows) and Program-3 (9 windows) were very 
similar. In some cases,, all the programs produced similar results (2 


8 Hz test) . 
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We can conclude that we will get narrow mainlobe widths (close 
to RSW) if we put the attenuation band close to the orign, but pay the 
price of low attenuation (-28 dB) . Conversely, good attenuation (-100 
dB) and broader mainlobes (twice that of RSW) are obtained for 
attenuation bands farther away from the origin. It also should be 
noted that opening up the bandwidth can have a stabilizing effect on 
the window's sidelobe behavior. A similar stability effect seems to be 
at work when we move the attenuation band away from the origin. 



21 


CHAPTER 4 

CONCLUSIONS AND FUTURE CONSIDERATIONS 
In this work we have reviewed the basic ideas involved with the 
Blackman- Tukey method of spectral estimation and, specifically, have 
investigated the problems associated with designing spectral lag 
windows. Three design ideas were implemented and tested. Program- 1 
used linear combinations of lag windows found in the literature. The 
combination coefficients were determined by solving a constrained 
least squares problem in which the objective function measured the 
spectral window response over a given frequency region. The 
constraints consisted of 1) nonnegativity on the combination 
coefficients and 2) that the combination coefficients sum to unity. 

For Program- 2, the constraints were changed so that the hybrid window 
itself was nonnegative. In Program- 3, the constraints of Program- 2 
were used with a different class of window functions (frequency 
variable cosinusoids) . 

Each of the three techniques was tested for performance over 
four different frequency bands, 1-2 Hz, 2 -4 Hz, 2 - 8 Hz, and 4 
8 Hz. In the case of Program- 3, the number of windows was varied to 
include the first three cosine windows, then the first six, and 
finally, the first nine. Based on these tests, we can tentatively 
conclude that: 

1) Relaxing the nonnegativity constraint on the a^ did not offer the 
advantages that were hoped for. The results for Program- 2 were at best 
comparable to Program- 1, and in some cases, not as good. 

2) The difference in behavior for Programs-1, 2 and 3, for attenuation 
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bands close to the origin and of a narrow bandwidth, is largely in how 
the sidelobes behave outside the attenuation band. The best results 
for this case were from Program-1 and Program-3 (3 windows). 

3) Using the cosine windows in Program- 3 to try to improve over the 
original idea of Program- 1 was a mixed success. Program- 3 seems mainly 
to offer advantages over Program- 1 in the case where the attenuation 
band is farther away from the origin (e.g., 4-8 Hz). For this case, 
Program- 3 achieved much better attenuation levels (at the cost of 
mainlobe width) provided we use enough windows. 

4) The level of attenuation improves and the main lobes get broader 
for all three programs as the attenuation band moves away from the 
origin. 

5) Increasing the attenuation bandwidth improves sidelobe behavior 
in general. We do not have to give up much in terms of attenuation 
levels and mainlobe widths. 

6) The number of cosine windows appropriate to use in Program- 3 varies 
as follows: for attenuation bands close to the origin, use fewer 
windows; for bands farther away, use more. Of course, the mainlobe 
widths will increase for the latter case. 

For future efforts along the line of research presented in this 
thesis, several things could be done. 

1) Implementation of the algorithm in which H is evaluated by 
numerical integration rather than by matrix multiplication (W W) may 
offer advantages. 

2) The three programs need to be studied more thoroughly in terms of 
their behavior with respect to various combinations of windows and 
spectral attenuation bands. 
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3) Modeling needs to be done to quantify the performance of the 
programs. For example, if we take two sinusoids of some relative 
magnitude at some frequency spacing, and truncate them to some 
specific length, we can then try to design a window that will resolve 
them. The relative strengths and frequency spacings could then be 
varied. 

4) In this study we ignored the sign of the hybrid spectral windows. 
The degree to which these windows have negative lobes and their 
characterization should be investigated. This should provide another 
criterion for judging window performance. It is desirable to have no 
negative sidelobes in the spectral window, but several typical window 
functions have some. The severity of this problem depends to some 
extent on the applications in mind. 
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TABLE 1 


Program- 1 and Program- 2 Lag Window 
Functions (defined on an interval [-1/2, 1/2] secs) 


WINDOW NO. 

NAME 

FORMULA w^r) 

1 

Rectangular 

1 

2 

Parzen- 2 

1 - 4r 2 

3 

Cosine- tip 

cos (nr ) 

4 

Bartlett-like 

1 + 2|r | 

5 

Hann 

.5 + .5 cos(2*r) 

6 

Hamming 

. 54 + .46 cos (2 nr) 

7 

Papoulis 1 - 2 1 r | 

|cos(2*r) + (1/n) | sin(2*r ) 

8 

Blackman ,42 + 

.5 cos(2*r) + .08 cos(47rr 

9 

Bartlett 

i 

ro 

10 

Sine -like 

sin(2*rr )/(2nr ) 

11 

Gaussian (a - 2.5) 

exp(-.5[2ar] 2 ) 
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TABLE 2 

Mean dB Levels (MDB) and Half Power Widths (HPW) 
for the 1 - 2 Hz Test 



MDB (db) 

HPW (Hz) 

PI 

-27 

0.48 

P2 

-28 

0.46 

P3 ; 3 

00 

CM 

i 

0.47 

P3;6 

-29 

0.46 

P3 ; 9 

-29 

0.46 
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TABLE 3 

Mean dB Levels (MDB) and Half Power Widths (HPW) 
for the 2 - 4 Hz Test 



MDB (db) 

HPW (Hz) 

PI 

-54 

0.66 

P2 

-58 

0.64 

P3 ; 3 

-55 

0.66 

P3 ; 6 

-60 

0.63 

P3 ; 9 

-64 

0.62 
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TABLE 4 

Mean dB Levels (MDB) and Half Power Widths (HPW) 
for the 2 - 8 Hz Test 



MDB (db) 

HPW (Hz) 

PI 

-56 

0.67 

P2 

-57 

0.67 

P3 ; 3 

-55 

0.68 

P3;6 

-57 

0.67 


P3 ; 9 


-57 


0.67 


28 


TABLE 5 

Mean dB Levels (MDB) and Half Power Widths (HPW) 
for the 4 - 8 Hz Test 



MDB (db) 

HPW (Hz) 

PI 

-83 

0.81 

P2 

-99 

0.85 

P3 ; 3 

-97 

0.85 

P3;6 

-110 

0.94 

P3 ; 9 

-117 

0.93 
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APPENDIX A 

In this appendix, we will show that the vector a that minimizes 
the unconstrained function a^Ha is a* — 0. The treatment is after 
Scales (1985). Consider a Taylor expansion for F(x) 

(1) F(x + Ax) - F(x) + Ax T g(x) + 1/2 Ax T G(x)Ax + . . . , 

where the gradient is g(x) - faF/dX}] 


LaF/3x n J 


and the Hessian is G(x) - \d 2 F/d* l dx 1 ... a 2 F/3x 1 ax n l 


_a 2 F/3x n ax 1 . . . 3 2 F/3x n 3x n J 


Now consider a quadratic form 


F(x) - 1/2 x T Ax + b T x + c. 


Then F(x + Ax) - 1/2 (x + Ax) T A(x + Ax) + b 1 (x + Ax) + c 


F(x) + Ax T (Ax + b) + 1/2 Ax a AAx . 


Comparing (2) with (1), we have that 


g (x) - Ax + b 


G(x) - A. 


The condition for x to be a stationary point is that 
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or 


g(x*) - 0 
Ax* - -b. 


Casting this result in the notation of our problem gives 

Ha* - -b. 

But, b - 0, 

so Ha* - 0. 

ic 

This linear system is consistent, with a solution a - 0. 
square, if it is nonsingular, the solution is unique, and 
solution. If H is positive definite it Is nonsingular. We 
require H to be positive definite although we only showed 
that It was nonnegative semidef inite . DQPROG will perturb 


Since H is 
0 is the 
actually 
in Chapter 2 
the 


diagonals so that H becomes positive definite If it is not. 
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APPENDIX B 

In this appendix we will first give some useful basic 
information about the specific quadratic programming (QP) subroutine 
we used to generate the thesis results, and then we will give a broad 
outline of the typical steps to solving a QP problem. 

In all three programs discussed in the body of this thesis the 
heart of each code was a subroutine from the IHSL MATH^-^- /LIBRARY 
called DQPROG. DQPROG is the double precision version of QPROG which 
is dicussed next. 

QPROG is based on M, J. D. Powell's implementation of the 
Goldfarb and Idnani (1983) quadratic programming algorithm for QP 
problems subject to general linear equality and inequality 
constraints. The matrix H (discussed in Chapter 2) is required to be 
positive definite. If it is not, then a positive definite perturbation 
is used in place of H. For more details, see Powell (1983a, 1983b). 

The following broad outline of quadratic programming is after 
Bronson (1982). The solution of many mathematical programming problems 
begins with the Kuhn-Tucker conditions. Consider the following general 
nonlinear programming problem 


given that 

X - [x lt . . 


maximize 

Z - f(x) 


subject to 

gi ^ 0, 

i - 1, m 

with 

X 

IV 

o 



To solve this program, rewrite the nonnegativity constraints as 


(iff! 1 
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S 0, 

and add slack variables to the left hand sides of the constraints. 
This converts the inequality constraints into equality constraints. 
Now form the Lagrangian function 

m m+n 

L - f(x) - S Ajjg^x) + (x n+1 ) Z ] - S Xil-y-i + C^+i) ], 
i-1 i-m+1 

where the A^ are Lagrange multipliers. Finally, solve the system 


dL/dXj - 0, 

j - 1. 

2n+m 

8L/dX i - 0, 

i - 1, 

m+n 

V 

(V 

o 

i - 1. 

m+n 


These last three expressions constitute the Kuhn-Tucker conditions. 
They are important because the set of their solutions contains the 
solution to our original nonlinear programming problem. 

The QP problem is a special case of the general nonlinear 
programming problem. Let's write the QP problem in the following form 
(note, seeking the maximum extremum is quite general in that the 
maxima of f(x) are the minima of -f(x)) 

maximize z - x^Cx + D^x, 

subject to Ax ^ B, 


with 


x £ 0. 


After applying the Kuhn-Tucker conditions to the QP problem, we find 
that the solution must satisfy the new matrix equation 
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where 


A'Y - B' , 


A I X 0 X 


-2C 0 


3 * I 2 


and 


S is a slack variable vector and U and V are Lagrange multiplier 
vectors. We also have the condition that 

Y' £ 0. 


Bronson (1982) discussed the use of the method of Frank and Wolfe to 
solve these equations. This method is rather complicated and beyond 
the scope of this thesis. However, in a few words, we can generally 
describe what takes place. A subproblera of our Kuhn-Tucker conditions 
is a modified linear programming (LP) problem with the inequality 
constraints of our QP problem. A well known method of solving LP 
problems is the Simplex method. A series of iterations that involves 
application of the Simplex method will give the desired solution if it 
exists . 


i 
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APPENDIX C 


PROGRAM PRGRAM1 


C Author: Kent Broadhead (12/89). 


C**SHt************************* 

INTEGER STRINGLEN .SIZE 

PARAMETER (IFFTMAX-2048 , IMAXCHNI^-50) 

INTEGER LDA , LDH , NCONMAX , NCON , NEQ , NVARMAX , NVAR 
PARAMETER (NCONMAX-12 , NEQ-1 , NVARMAX-11 
1 , LDA-NCONMAX , LDH-NVARMAX) 

CHARACTER* 20 , FILEOUT , DIGNSTC , TIME 
CHARACTER QUESTN*1 
DIMENSION IARRAY(ll) 

REAL* 8 AME (IMAXCHNL, NVARMAX) , 

1 AMET (NVARMAX , IMAXCHNL) 

DIMENSION SCALE (NVARMAX) , SUMARRY(IFFTMAX) 

1 , SUMDB ( I FFTMAX) 

REAL* 8 A(LDA , NVARMAX) ,ALAMD A (NVARMAX) 

I , B (NCONMAX) , DIAG , G (NVARMAX) , 

1 H(LDH,LDH) , SOL(NVARMAX) 

COMPLEX CW( I FFTMAX) 

INTEGER BEGCHNL , ENDCHNL 
INTEGER IACT (NVARMAX) , K , NACT , NOUT 
EXTERNAL DQPROG , UMACH 

C*****************'******'************ 


PRINT *, 'Enter the output file name:' 

READ (*,' (A20)') FILEOUT 

OPEN ( 5, FILE-FILEOUT, STATUS- 'NEW' ) 


C - Window only option 

PRINT *,'Do you just want a window plot? 

1 (Y or RETURN) : ' 

READ ( * , ' ( Al ) ' ) , QUESTN 

IF (QUESTN .EQ. 'Y' .OR. QUESTN .EQ.'y') THEN 
CALL WNDONLY(CW , SUMDB , ILAG1 , SUMARRY) 

GO TO 1050 
ENDIF 

C 

C 

C Get the main input params. No values for the arguments 
C are sent. Both returned by GETLIST. 

CALL GETLIST (I ARRAY, NVAR) 

WRITE (*, '(IX, 12, 4X, 12)'), ( IAKRAY ( I ) , I , 1-1 , NVAR) 


PRINT *, 'Enter no. of windows to use (2-11):' 
ACCEPT *, NVAR 
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C SCALE(i) sums to unity constraint 
DO 302 J-l.NVAR 
302 A(1,J)-1.0 

C Open diagnostics file; 

SIZE-STRINGLEN ( FILEOUT , LEN ( FILEOUT) ) 
DIGNSTC - FILE0UT(1:SIZE)//' .DIAG' 

OPEN (4 , FILE-DIGNSTC , STATUS- 'NEW' ) 
WRITE(4,*) , 'Diagnostics File' 
WRITE(4,*) ' ' 

WRITE (4, *) , 'VAX filename - PRGRM1 . FOR' 
WRITE(4 , *) ' ' 

WRITE (4 ,444) , FILEOUT 
444 FORMAT ( ' Data filename - ' ,A20) 
WRITE(4 ,*) ' ' 

C Windows range from - . 5 to .5 
C Get 0 to . 5 first, then use symm to pad 
C out neg side (in FFT format). 


C Initialize 

ICNT - 1 

PRINT *, 'No. smples for lag wndw (desgn) ' 
READ(*,*) , ilagl 

PRINT *, 'No. smples for spec wndw (desgn)' 
READ(* , *) , ifftnml 

XLAG1 - ILAGl 
XFFTNM1 - IFFTNM1 
DELTAF1 - XLAG 1 /XFFTNM1 

PRINT *, 'No. smples for lag wndw (dsply) ' 
READ(*,*) , ilag2 

PRINT *, 'No. smples for spec wndw (dsply)' 
READ(* ,*) , ifftnm2 

XLAG 2 - ILAG2 
XFFTNM2 - IFFTNM2 

DELTAF2 - XLAG2/XFFTNM2 

ihalfl - ilagl/2 + 1 
ihalf2 - ilag2/2 + 1 
NYQstl - ifftnml/2 + 1 
NYQST2 - iff tnm2/2 + 1 

WRITE(4, 1010) .ilagl 

1010 FORMAT ( IX ,' No . smples for lag wndw (desgn)', 14) 
WRITE(4, 1011) .ifftnml 

1011 FORMAT (IX, 'No. smples for spec wndw (desgn) ',14) 


WRITE(4 , 1020) , ilag2 

1020 FORMAT ( IX ,' No . smples for lag wndw (dsply)',I4) 
WRITE(4 , 1021) , ifftnm2 

1021 FORMAT ( IX No . smples for spec wndw (dsply)',I4) 
WRITE(4 , *) ' ' 

PRINT *, 'Enter the begining and ending freqs:' 
c PRINT *, 'Note: max no. of channels 

C 1 presently is 50.' 

ACCEPT *,F1,F2 

C Cnvert FI and F2 to fft points. 

CALL CNVRT ( BEGCHNL , ENDCHNL , FI , F2 , DELTAF1 ) 

WRITE(4 , 1000) , FI , F2 

1000 FORMAT ( IX, 'Frq atten. band \F10.3,' to 'F10.3) 
WRITE(4 , *) ' ' 

WRITE(4 , 1023) , BEGCHNL , ENDCHNL 
1023 FORMAT ( IX, 'Chnl band ',i4,' to 'i4) 

WRITE(4 , *) ' ' 

NCON - NVAR + 1 

C Set up constraints; 

DO 305 ICLR— 1 , NVARMAX 

305 G(ICLR)-0.0D0 

B(l) - 1.0D0 
DO 299 1-2, NCON 
299 B(I) - 0.0D0 

DO 301 I-l.NCON 

DO 301 J-l.NVAR 
301 A(I , J)— 0 . 0D0 

DO 307 J-l.NVAR 

307 A(1 , J)— 1 . 0D0 

DO 303 I-l.NVAR 

303 A(I+1 , I)— 1 . 0D0 

C Begin window calcs : 
xlagl - ilagl 
DO 20 I OUTER-1 , NVAR 

CALL CCLEAR(CW, ifftnml) 

C Compute a window. 

DO 10 1-1 , ihalfl 
XI - I 

ARG - (XI- 1 . ) /xlagl 

10 CW(I)-WIND0W(ARG , IARRAY(ICNT) ) 

C Fourier transform the window. 

CALL FFT (CW, ifftnml) 

C Load the window transform as a column of AME 

DO 40 IL-BEGCHNL, ENDCHNL 
INDEX - IL-BEGCHNL+1 
40 AME ( INDEX, ICNT) - CW(IL) 

ICNT - ICNT + 1 
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CONTINUE 


C Define M. 

M - ENDCHNL- BEGCHNL+1 

C Take transpose. 

DO 12 K-l.M 

DO 50 L-l.NVAR 

50 AMET(L,K) - AME(K,L) 

12 CONTINUE 

C Multiply transpose by orig matrix and store in ATA. 
DO 100 I-l.NVAR 
DO 100 J-1,NVAR 
SUM-0 . 

DO 102 K-l.M 

102 SUM - SUM + AMET ( I , K) *AME (K , J ) 

100 H(I,J)- 2.*SUM 

C DO 13 K-l.NVAR 

C WRITE (6 , # (IX, 8F10 . 5) ' ) (H(K,L) , L-1,NVAR) 

C13 CONTINUE 


C Call DQPROG to solve for the hybrid window scale factors; 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IMSL Name: 
Computer : 
Revised: 
Purpose : 

Usage : 


QPROG /DQPROG (Single/Double precision version) 

VAX/SINGLE 

October 15, 1985 

Solve a quadratic programming problem subject 
to linear equality/inequality constraints . 

CALL QPROG (NVAR, NCON, NEQ, A, LDA, B, G, H, 
LDH.DIAG, SOL, NACT, I ACT, ALAMDA) 


Arguments : 

NVAR - The number of variables. (Input) 

NCON - The number of linear constraints. (Input) 

NEQ - The number of linear equality constraints. (Input) 

A - Real NCON by NVAR matrix. (Input) 

The matrix contains the equality contraints in the 
first NEQ rows, followed by the inequality 
constraints . 

LDA - Leading dimension of A exactly as specified in the 
dimension statement of the calling program. (Input) 
B - Real vector of length NCON containing right-hand 

sides of the linear constraints. (Input) 

G - Real vector of length NVAR containing the 

coefficients of the linear term of the objective 
function. (Input) 

H - Real NVAR by NVAR matrix containing the Hessian 

matrix of the objective function. (Input) 

H should be symmetric positive definite, if H 
is not positive definite, the algorithm attempts 
to solve the QP problem with H replaced by a 
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C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


LDH 

DIAG 

SOL 

NACT 

IACT 

ALAMDA 


H + DIAG* I such that H + DIAG*I is positive 
definite - See Remark 3. 

Leading dimension of H exactly as specified In the 
dimension statement of the calling program. (Input) 
Real scalar equal to the multiple of the 
identity matrix added to H to give a positive 
definite matrix. (Output) 

Real vector of length NVAR containing solution. 
(Output) 

Final number of active constraints. (Output) 

Integer vector of length NVAR containing the indices 
of the final active constraints in the first 
NACT positions. (Output) 

Real vector of length NVAR containing the 
Lagrange multiplier estimates of the final 
active constraints in the first NACT positions. 
(Output) 


CALL DQPROG (NVAR , NCON , NEQ , A , LDA , B , G , H , LDH , DIAG , SOL , 
1 NACT, IACT, ALAMDA) 


C CALL UMACH ( 2 , NOUT ) 

PRINT * , ' 1 

WRITE(*,*) ( SOL(K) , K-l , NVAR) 

C WRITE(*, 99999) (SOL(K) , K-l , NVAR) 

C99999 FORMAT (IX, 8F10 . 5) 

SUM-0. 

DO 432 1-1 , NVAR 

SCALE (I) - SOL(I) 

432 SUM - SUM + SOL(I) 

WRITE (*, ' (1X,4HSUM— , F10 . 5) ' ) SUM 


C Put some Info in the diag. file; 

WRITE (4, *) , 'Optimum Window Wts' 

WRITE(4 , *) (SCALE(I) , 1-1 ,NVAR) 

WRITE(4 , *) ' ' 

WRITE (4,*) 'The program uses the following percentages 
1 for the given windows ; ' 

DO 452 I PRC NT -1 , NVAR 

PERCNT - SCALE(IPRCNT)*100 . 

WRITE(4,453) , IARRAY(IPRCNT) ,PERCNT 
453 FORMAT ( 9 WINDOW NO . ',13 , 1 - ',F10.2,' PERCNT') 

452 CONTINUE 

WRITE(4 , *) ' ' 

C - 

C Now that we have the lin. comb, wts, we need to get 
C the spectrum of the hybrid window. To do this, let's 
C recompute the window functions, scale them, 

C stack them, and FFT. 


C Initialize 
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ICNT - 1 

CALL CLEAR (SUMARRY, ifftnm2) 

C Begin window calcs: 
xlag2 - ilag2 
DO 22 I OUTER-1 , NVAR 

C Compute, scale, and stack a window. 

DO 11 1-1 , ihalf 2 
XI - I 

ARG - (XI -1 . )/xlag2 

11 SUMARRY(I)— SUMARRY(I)+SCALE(ICNT)* 

1 WINDOW (ARG, ICNT) 

ICNT - ICNT + 1 
22 CONTINUE 

C Put hybrid window (time domain) out to a plot file. 
SIZE— STRINGLEN(FILEOUT , LEN(FILEOUT) ) 

TIME - FILEOUT(l: SIZE)// 'T. DAT’ 

OPEN ( 3 , FILE-TIME , STATUS- ' NEW ' ) 

THALF - 0.5 

XHALF2 - IHALF2 

DELTAT - THALF/ (XHALF2 - 1.) 

SUM - 0.0 
T - 0.0 

WRITE (3 ,*) T , SUMARRY(l) 

DO 405 J 5-2, ihalf 2 
SUM - SUM + DELTAT 
T - SUM 

WRITE (3 ,*) T, SUMARRY (J5) 

405 CONTINUE 


C Fourier transform the window. 

CALL CCLEAR(CW, ifftnm2) 

DO 23 IMOV-l , ihalf2 

23 CW(IMOV) -SUMARRY (IMOV) 

CALL FFT(CW, If f tnm2) 

CALL CLEAR (SUMARRY, iff tnm2) 

DO 24 IMOV-l, NYQST2 

24 SUMARRY ( IMOV )-CW( IMOV) 

C Output Stage: 

C Basic set up 

XMAX - SUMARRY ( 1 ) 

DO 400 IDB-1 .NYQST2 

ARG - ABS ( ( SUMARRY( IDB) /XMAX) ) 

IF (ARG. EQ. 0.0) THEN 

PRINT *,'A zero value was encountered in 
1 the hybrid spectrum window.' 

PRINT *,'The dB value has been set to -100 
1 0000000. I recommend that' 

PRINT *,'you go with the clipping preset on 



1 the next question. ' 

SUMDB(IDB) - -1000000000. 

ELSE 

SUMDB(IDB) - 20 . *AL0G10(ARG) 

ENDIF 
400 CONTINUE 



C Find mean of attenuated zone; 

C First, get chnnel pts. for given freq band; 

CALL CNVRT(BEGCHNL, ENDCHNL, FI , F2 , DELTAF2 ) 

XSUM - 0. 

DO 402 IMEAN-BEGCHNL , ENDCHNL 
402 XSUM - XSUM + SUMDB(IMEAN) 

XNUMBR - ENDCHNL - BEGCHNL + 1 
XMEANDB - XSUM / XNUMBR 
WRITE (4, 406) XMEANDB 
406 FORMAT (IX, 'XMEANDB -' ,G20. 10) 

WRITE(4 ,*) ' ' 

C Find Half Power Width: 

CALL GETHFPW ( SUMDB , HFPW , DELTAF2 ) 

WRITE (4, 873) ,HFPW 

873 FORMAT (IX, 'Half Power Freq - ' , F6 . 2 , IX, 'Hz ' ) 
WRITE(4 , *) ' ' 

C Plot file details 

PRINT *, 'Clipping? (N or Return):' 

READ(* , ' (Al)') .QUESTN 

IF (QUESTN .EQ. 'N' .OR. QUESTN .EQ.'n') THEN 
C Output data as is . 

DO 404 J-1.NYQST2 
FREQ - (J-1)*DELTAF2 

WRITE (5 ,*) FREQ, SUMDB (J) 

404 CONTINUE 

ELSE 

C Output data with clipping 

CLIP - -160.0 

PRINT *, 'Enter Y (change clipping value) 
lor Return (for -160 preset):' 

READ (*,'(A1)'), QUESTN 

IF (QUESTN .EQ. 'Y' .OR. QUESTN .EQ.'y') THEN 
PRINT *,' Enter new clipping value:' 
READ(*,*) .CLIP 
ENDIF 

DO 401 J-1.NYQST2 

FREQ - (J-1)*DELTAF2 
IF(SUMDB(J) .GE.CLIP) THEN 
WRITE (5,*) FREQ.SUMDB(J) 

ELSE 

WRITE (5,*) FREQ, CLIP 


o o 


401 


END IF 
CONTINUE 


ENDIF 
C 


C Put some info in the diag. file; 

WRITE (4, *), 'The first few values' 
WRITE (4 , *) , 'of the plot file:' 

DO 403 J2-1.10 
XJ2 - J2-1 

WRITE(4 , *)XJ2 , SUMDB(J2) 

403 CONTINUE 

C END OF MAIN PROGRAM 
1050 CONTINUE 
END 


BEGIN SUBROUTINES 

SUBROUTINE GETHFPW(X , F2 , DELTAF) 
DIMENSION X(l) 


C Find dB pt. just past -3dB. 

DO 10 1-1,10000 

IF(X(I) .LE. -3.) THEN 
IF(X(I) .EQ. -3.) THEN 

12 - I 

CALL CNVRT2( I, F2, DELTAF) 

GO TO 20 
ELSE 

II - 1-1 

13 - I 

CALL CNVRT2 ( II , FI , DELTAF) 

CALL CNVRT2(I3 , F3 .DELTAF) 

CALL INTERP(F1 ,F2,F3,X(I1) ,X(I3)) 
GO TO 20 
ENDIF 
ENDIF 

10 CONTINUE 

20 CONTINUE 

RETURN 

END 

SUBROUTINE CNVRT2( I, F, DELTAF) 


F - ( I - 1 ) +DELTAF 


RETURN 

END 

SUBROUTINE INTERP (FI , F2 , F3 , DB1 ,DB3) 
DB2 - -3. 


C Convert dB to ratio; 

R1 -( 10.**(DBl/20. )) 
R2 -( 10.**(DB2/20. )) 
R3 -( 10 . **(DB3/20 . ) ) 
C Linearly interpolate; 


F2 - FI + (F3-F1)*(R2-R1)/(R3-R1) 

RETURN 

END 

SUBROUTINE CNVRT ( BEGCHNL , ENDCHNL , FI , F2 , DELTAF) 
INTEGER BEGCHNL, ENDCHNL 

BEGCHNL - (Fl/DELTAF) +1.5 
ENDCHNL - (F2/DELTAF) +1.5 

I WIDTH - ENDCHNL - BEGCHNL + 1 
IF (I WIDTH .GT. 50) THEN 
PRINT *, 'Channel width .GT. 50' 

STOP 
END IF 
RETURN 
END 

SUBROUTINE FFT(CW.LX) 

COMPLEX CW(1) 

CALL FFTWRAP ( CW , LX) 

CALL FORK(LX,CW, 1.0) 

RETURN 

END 

SUBROUTINE FFTWRAP (CW, LX) 

COMPLEX CW(1) 

NYQUI ST-LX/2 + 1 
I END-LX - NYQUI ST 
DO 10 1-1 , I END 

10 CW( NYQUI ST+I) - CW(NYQUIST-I) 

RETURN 

END 

FUNCTION WINDOW(ARG.ICNT) 

PI - 3.14159 
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IF(ICNT.EQ.l) THEN 
GO TO 1 

ELSEIF ( ICNT . EQ . 2 ) THEN 
GO TO 2 

ELSEIF(ICNT . EQ . 3) THEN 
GO TO 3 

ELSEIF(ICNT . EQ . 4) THEN 
GO TO 4 

ELSEIF(ICNT . EQ . 5) THEN 
GO TO 5 

ELSEIF (ICNT . EQ . 6 ) THEN 
GO TO 6 

ELSEIF (ICNT . EQ . 7) THEN 
GO TO 7 

ELSEIF(ICNT . EQ . 8) THEN 
GO TO 8 

ELSEIF(ICNT.EQ.9) THEN 
GO TO 9 

ELSEIF(ICNT.EQ.IO) THEN 
GO TO 10 

ELSEIF(ICNT.EQ.ll) THEN 
GO TO 11 
ELSE 

GO TO 100 

END IF 

CONTINUE 

Rectangular 

WINDOW - 1. 

GO TO 100 

CONTINUE 

Parzen-2 

WINDOW - 1. - 4. * (ARG**2) 

GO TO 100 

CONTINUE 

Cosine-Tip 

WINDOW - COS (ARG*PI) 

GO TO 100 

CONTINUE 

Bartlett -1 ike , with sign change 
WINDOW - 1. + 2.*ABS(ARG) 

GO TO 100 

5 CONTINUE 
C Hann 

WINDOW - .5 + . 5*COS(2.*ARG*PI) 
GO TO 100 

6 CONTINUE 
C Hamming 
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WINDOW - .54 + ,46*COS(2.*ARG*PI) 

GO TO 100 

7 CONTINUE 
C Papoulis 

WINDOW - 1. - 2 .*ABS(ARG)*COS(2.*ARG*PI) 

1 + ( 1 . /PI ) *ABS ( SIN( 2 . *PI*ARG) ) 

GO TO 100 

8 CONTINUE 
C Blackman 

WINDOW - .42 +.5*COS(2.*ARG*PI)+.08*COS(4.*ARG*PI) 

GO TO 100 

9 CONTINUE 
C Triangle 

WINDOW - 1. - 2 . *ABS (ARG) 

GO TO 100 

10 CONTINUE 

IF (ARG .EQ. 0.) THEN 
WINDOW - 1. 

ELSE 

WINDOW - SIN(2.*PI*ARG)/(2 .*PI*ARG) 

ENDIF 
GO TO 100 

11 CONTINUE 
C Gaussian 

ALPHA - 2.5 

ARG2 - - .5*((2.*ALPHA*ARG)**2) 

WINDOW - EXP (ARG 2) 

GO TO 100 

100 CONTINUE 

RETURN 

END 

SUBROUTINE FORK(LX,CX, SIGNI) 

C From Claerbout's FGDP (first ed.), pg 12. Typed in by MKB. 

COMPLEX CX(LX) , CARG , CEXP , CW, CTEMP 
J-l 

SC-SQRT ( 1 . /LX) 

DO 30 I-l.LX 
IF(I.GT.J) GO TO 10 
CTEMP— CX ( J ) * S C 
CX(J)-CX(I)*SC 
CX(I)— CTEMP 
10 M— LX/2 

20 IF(J.LE.M) GO TO 30 

J-J-M 


IF(M.GE.l) GO TO 20 

J-J+M 

L-l 

ISTEP-2*L 
DO 50 M-1,L 

CARG— (0 . , 1 . )*(3 . 14159265*SIGNI*(M-1) )/L 
CW-CEXP ( CARG ) 

DO 50 I— M, LX, I STEP 
CTEMP-CW*CX(I+L) 

CX ( I+L) -CX ( I ) - CTEMP 

CX(I)-CX(I)+CTEMP 

L-ISTEP 

IF(L.LT.LX) GO TO 40 

RETURN 

END 

SUBROUTINE CCLEAR(X.LEN) 

COMPLEX X(l) 

DO 10 I-l.LEN 

X(I) - (0..0.) 

RETURN 

END 

SUBROUTINE CLEAR (X.LEN) 

DIMENSION X(l) 

DO 10 I-l.LEN 
X(I) - 0. 

RETURN 

END 


integer function stringlen(string, length) 
integer length 
character*(*) string 
logical done 
integer i 
i-length 
done— . false . 
do while (.not. done) 

done— (string( i : i) .ne . ' ' ) .or . (i. le . 1) 

i-i-1 

end do 


check for "empty" string; i.e., having only blanks 
if (i . eq . 1 .and. string(l:l) .eq. ' ') i-0 

stringlen-i 

if(i.gt.O) stringlen-i+1 

return 

end 


SUBROUTINE WNDONLY ( CW , SUMDB , ILAG1 , SUMARRY) 


o o 
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COMPLEX CW(1) 

DIMENSION WIND (1024) , SUMDB(l) , SUMARRY( 1) 

CHARACTER 0PTI0N*4 
CHARACTER QUESTN*1 
PRINT * , ' Enter window no . : ' 

READ(* , *) , IWINDNO 
C Calculate window. 

PRINT *, 'No. smples for lag wndw (desgn) ' 

READ(*,*) , ilagl 
CALL CALC (WIND , IWINDNO , ILAG1) 

C Prepare output. 

PRINT *, 'Enter FUNC or FREQ:' 

READ(* , 10) .OPTION 
10 FORMAT (A4) 

IF(OPTION .EQ. 'FUNC') THEN 

Output function domain window values for plotting. 
CALL FUNCOUT (WIND , ILAGl) 

ELSE 

Compute FFT and output freq. domain window values 
for plotting. 

PRINT *, 'No. smples for spec wndw (desgn)' 

READ(* ,*) , ifftnml 

CALL FREQOUT (WIND , CW , SUMDB , SUMARRY 
1 , ILAGl , IFFTNM1) 

ENDIF 

RETURN 
END 

SUBROUTINE CALC (WIND , IWINDNO , ILAG1) 

DIMENSION WIND(l) 
xlagl - ilagl 
ihalfl - ilagl/2 + 1 
C Compute a window. 

DO 10 1-1, ihalfl 
XI - I 

ARG - (XI -1.) /xlagl 
10 WIND(I)-WINDOW(ARG, IWINDNO) 


RETURN 

END 

SUBROUTINE FUNCOUT (WIND , ILAG ) 
DIMENSION WIND(l) 

IHALF2 - I LAG/2 + 1 

THALF - 0.5 

XHALF2 - IHALF2 

DELTAT - THALF/ (XHALF2 - 1.) 

SUM - 0.0 

T - 0.0 

WRITE (5 ,*) T.WIND(l) 

DO 405 J5-2 , ihalf2 
SUM - SUM + DELTAT 




aiiml mj |, | 


129 


T - SUM 

WRITE (5 ,*) T,WIND(J5) 

405 CONTINUE 

RETURN 

END 

SUBROUTINE FREQOUT (WIND , CW , SUMDB , SUMARRY 
1 , ILA.G2 , IFFTNM2 ) 

DIMENSION WIND(l) , SUMARRY(l) , SUMDB(I) 
COMPLEX CW(1) 

CHARACTER QUESTN*1 
NYQST2 - IFFTNM2/2 + 1 
XLAG2 - I LAG 2 
XFFTNM2 - IFFTNM2 
DELTAF2 - XLAG2/XFFTNM2 
IHALF2 - ILAG2/2 + 1 

CALL CCLEAR(CW,ifftnm2) 

DO 23 IMOV-l , ihalf2 

23 CW ( IMOV) —WIND ( IMOV) 

CALL FFT(CW, Ifftnm2) 

CALL CLEAR( SUMARRY, if ftnm2) 

DO 24 IMOV-1, NYQST2 

24 SUMARRY ( IMOV )-CW( IMOV) 


C Output Stage: 

C Basic set up - 

XMAX - SUMARRY ( 1 ) 

DO 400 IDB-1 , NYQST2 

ARG - ABS( ( SUMARRY (IDB) /XMAX)) 

IF (ARG. EQ. 0.0) THEN 

PRINT * , ' A zero value was encountered in 
1 the hybrid spectrum window.' 

PRINT *,'The dB value has been set to -100 
1 0000000. I recommend that' 

PRINT *,'you go with the clipping preset on 

1 the next question.' 

SUMDB (IDB) - -1000000000. 

ELSE 

SUMDB(IDB) - 20.*ALOG10(ARG) 

ENDIF 
400 CONTINUE 




C Find Half Power Width: 

CALL GETHFPW ( SUMDB , HFPW , DELTAF2 ) 
WRITE(* , 873) ,HFPW 
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873 

C--- 

C 

404 

C 


401 


FORMAT (IX, 'Half Power Freq - ' , F6 . 2 , IX, 'Hz ' ) 


Plot file details - 

PRINT *, 'Clipping? (N or Return):' 

READ(* , ' (Al)') .QUESTN 

IF (QUESTN .EQ. 'N' .OR. QUESTN .EQ.'n') THEN 
Output data as is . 

DO 404 J-l ,NYQST2 
FREQ - (J-1)*DELTAF2 

WRITE (5,*) FREQ, SUMDB(J) 

CONTINUE 

ELSE 

Output data with clipping 
CLIP - -160.0 

PRINT *, 'Enter Y (change clipping value) 
lor Return (for -160 preset):' 

READ(* , ' (Al) ' ) .QUESTN 

IF (QUESTN .EQ. 'Y' .OR. QUESTN .EQ.'y') THEN 
PRINT Enter new clipping value:' 

READ (*,*) .CLIP 
ENDIF 

DO 401 J-l , NYQST2 

FREQ - (J -1)*DELTAF2 
IF(SUMDB(J) .GE.CLIP) THEN 
WRITE (5,*) FREQ, SUMDB(J) 

ELSE 

WRITE(5 ,*) FREQ, CLIP 
ENDIF 
CONTINUE 
ENDIF 

RETURN 

END 

SUBROUTINE GETLIST(IARRAY.NUM) 

DIMENSION I ARRAY ( 1 ) 

INTEGER NUM 

PRINT *,' ' 

PRINT *, 'Welcome to PRGRM1 , an alternative window 
1 design program' 

PRINT *,' ' 

PRINT *, 'Later below, you will be asked for a list 
1 of numbers corresponding to' 

PRINT *,'the windows you want to use. The window options 

1 are given here with' 

PRINT *, 'their respective identifying numbers.' 

PRINT *,' ' 

PRINT *, 'Window options:' 

PRINT *, 'Window No. 1 - Rectangular' 

PRINT *, 'Window No. 2 - Parzen-2' 
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PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
PRINT *, 
ACCEPT * 


'Window No. 3 - 
'Window No. 4 - 
'Window No. 5 - 
'Window No. 6 - 
'Window No. 7 - 
'Window No . 8 - 
'Window No. 9 - 
'Window No. 10 
'Window No. 11 


Cosine-Tip' 

Bartlett- like, with sign change 
Hann' 

Hamming' 

Papoulis' 

Blackman' 

Bartlett (Triangle)' 

Sine like' 

Gaussian' 


'Enter the total number of windows you are' 
'going to use (between 2 and 11 inclusive):' 

, NUM 


C Do some editing. 

IF (NUM .GE. 2 


.AND. NUM . LE. 11) THEN 


ELSE 
PRINT *, 
PRINT *, 
STOP 
END IF 


'The number of windows you can use must' 
'be between 2 and 11 (inclusive). 


PRINT * ' ' 

PRINT *, 'Enter window no.s here (list of integers, each on 
PRINT *, 'a separate line):' 

READ(* , 100) , (IARRAY(I),I-1,NUM) 

100 FORMAT (12) 

C Do some editing. 

DO 30 I-l.NUM 

IF (IARRAY(I) .GE. 1 .AND. IARRAY(I) .LE. 11) THEN 
ELSE 

PRINT * , ' One or more of the window no.s you gave is outside' 
PRINT *, 'the acceptable range (1-11).' 

STOP 

ENDIF 

30 CONTINUE 


END 
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